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COMBINATION VERSUS SEQUENTIAL MONOTHERAPY IN CHRONIC 
HBV INFECTION: A MATHEMATICAL APPROACH 

DANIELA BERTACCHI, FABIO ZUCCA, SERGIO FORESTI, DAVIDE MANGIONI, AND ANDREA GORI 



Abstract. Sequential monotherapy is the most widely used therapeutic approach in the treatment 
of HBV chronic infection. Unfortunately, under therapy, in some patients the hepatitis virus mutates 
and gives rise to variants which are drug resistant. We conjecture that combination therapy is able 
to delay drug resistance for a longer time than sequential monotherapy. To study the action of 
these two therapeutic approaches in the event of unknown mutations and to explain the emergence 
of drug resistance, we propose a stochastic model for the infection within a patient which is treated 
with two drugs, either sequentially or contemporaneously, and develops a two-step mutation which 
is resistant to both drugs. We study the deterministic approximation of our stochastic model and 
give a biological interpretation of its asymptotic behaviour. We compare the time when this new 
strain first reaches detectability in the serum viral load. Our results show that the best choice is 
to start an early combination therapy, which allows to stay drug-resistance free for a longer time. 
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Hepatitis B virus is the cause of one of the most common infections in the world, with around 400 
million people infected according to the WHO assessment [30] • Its consequences, such as the liver 
cirrhosis and the hepatocellular carcinoma, affecting approximately 25% of patients with chronic 
HBV infection, are important causes of morbidity and mortality worldwide [30 1 126 |. ITT, I34j . 

The progression of the liver disease relies on environmental factors (co-infection with other viruses 
such as HCV or HIV, alcohol abuse), host factors (age, immune status), and viral factors. As to the 
latter, a great importance is assumed by the viral load and the presence of HBV genomic mutations, 
that can significantly affect the response to the treatment [33 J . Antiviral therapy has demonstrated 
to have an important role in reducing both the markers of liver disease progression (with new 
generation nucleoside analogues, viral suppression, i.e. the undetect ability of HBV DNA, can be 
reached in about 95% of all the cases [16]) and, above all, the evolution to cirrhosis and hepatocel- 
lular carcinoma |34j . It has been shown that high serum levels of HBV-DNA (> 10 4 copies/ml) are 
predictive of an increased risk of hepatocellular carcinoma, regardless of the transaminase level and 
the presence of liver cirrhosis [11] : antiviral therapy should therefore be considered a preventive 
and anti-cancer therapy [34J. 

Over recent years several antiviral drugs specific for the virus enzyme DNA polymerase have 
entered in the clinical practice, but so far they are mainly used in a sequential monotherapy (that 
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is, only one drug is employed at a time and if an increase in viral load is observed, then the 
drug is substituted with another one). Moreover, interferon monotherapy is still a mainstay for 
mild to moderate chronic infections [20J. Some unsolved problems remain today: the clearance, 
or viral eradication (certified by the disappearance of HBsAg or cccDNA), which would allow 
the discontinuation of the therapy, is obtained in only 10% of treated patients [27 1; to improve 
this result, different steps of HBV replicative cycle are being considered as possible target for 
the development of new drugs [57]. On the other hand, the long term treatment with sequential 
monotherapy [2U\ [TU\ I46| . is likely to select drug resistant strains, because of the onset of viral 
breakthrough and mutant escape (that is the emergence of variant viruses) under the selective 
pressure of medications [23J IS] ■ 

As stated in the international guidelines, the combination therapy (that is, therapy with a 
cocktail of drugs) has a well-defined role in the treatment of chronic hepatitis B, especially as 
salvage therapy and first line application in selected categories of patients, such as those with 
decompensated cirrhosis, HIV co-infection, pre-existing mutations and post-liver transplantation 
[2DJ QUI S3 EE] • The guidelines, however, are based only on published results, and to date, only a few 
studies have attempted to compare the effectiveness of the two different approaches (combination 
versus monotherapy) in treatment-naive patients. These studies did show promising results in favor 
of combination therapy ( [H1I511 I2T1 [53] ) . 

We conjecture that combination therapy, as a therapy for chronic HBV infection, is preferable 
to sequential monotherapy. Our conjecture is based on the following observation: despite being a 
stable DNA virus, hepatitis B virus replicates in the host cell through a RNA phase |37j . that makes 
it similar to the Human Immunodeficiency Virus (HIV). The natural history of HIV infection has 
been radically transformed by the introduction of combination antiviral therapy (Highly Active 
Antiretroviral Therapy, HAART) [9], therefore we believe that this kind of therapy should be 
beneficial also for HBV treatment. In particular we infer that under combination therapy the 
appearance of drug resistance is less likely, or at least, occurs after a longer time. We believe it 
is important to compare the outcome of combination therapy and of sequential monotherapy in 
the treatment of chronic HBV infection using a mathematical model, which takes into account the 
emergence of mutants of the virus. 

During the last two decades and more, mathematical models have been proposed to describe 
viral infections (such as HBC, HCV, HIV) both from a single patient's point of view, which is our 
point of view in this paper, and from an epidemics point of view. The standard mathematical 
model for chronic HBV infection within a single patient is given by a three ordinary differential 
equations system (in short, ODEs system): 

' f = A - S'T - aVT; 

Y = aVT-5Y; (1.1) 

V = P Y - cV; 



where T denotes target (uninfected) hepatocytes which are produced at rate A, die at rate 5'T and 
are infected at a rate aVT; Y denotes infected hepatocytes which die at rate 5Y (5 possibly larger 
than 5') and V denotes free virus. Virions are produced from infected cells at rate pY and are 
cleared from bloodstream at rate cV. The system has no exact solutions but approximate solutions 
may be obtained under different hypotheses. 

If the patient undergoes treatment, the drug acts against reproductions and/or infections. In 
the model this is explained by introducing the efficacy parameters e € [0, 1] ([331 [52]) and r\ G [0, 1] 
([3SJ|32]) and replacing the previous system by 

't = X - 5'T - (1 - rj)aVT; 

Y = (l-rj)aVT-5Y; (1.2) 

V = {l-s)pY -cV. 

The parameters e and rj respectively indicate the ability of treatment to block viral production or 
infections (the extremal values and 1 refer to no ability at all and complete ability, respectively). 
Of course the behaviour under no treatment is retrieved setting rj = e = 0. 

This model has been used in several papers to fit experimental data under treatment with 
different drugs, leading to estimates for some of the parameters involved, mainly for e, c and 8 
( [381 EH [561 EH [321 E21 H2| ) . It has also been modified to model HIV infection ([13]) or acute HBV 
infection ([14]). 

To our knowledge, none of the existing models gives an explanation of the emergence of drug- 
resistant strains of the virus, which are usually not observed without therapy. We want to modify 
the previous model, according to some biological observations. Within a single infected person, 
genetically different viral variants can coexist: that is because of the same nature of HBV, which 
reproduces itself through a pre-genomic RNA phase, accounting for the presence of heterogeneous 
strains in the viral progeny |371 144] . The mutations are pre-adaptive, linked to random errors of the 
polymerase enzyme, especially in the RNA phase of the replicative process. Hence, the antiviral 
therapy does not determine any resistant strains, but simply selects them |37 1 144 } [35 | I40| . Moreover, 
in treatment-naive patients a minority of drug-resistant species can also be found. These variants 
could be the result of the persistence of drug-resistant viruses transmitted by a treated individual, 
or may have been generated ex novo during the course of the original infection [44] . 

In order to keep the model simple and to focus on the first mutations which appear under 
sequential monotherapy, we think of the example of a treatment-naive patient which undergoes 
sequential monotherapy, with two drugs at hand (drug A and drug B) and who will develop drug- 
resistance to both. Before therapy, he/she has a high serum viral load of HBV virions, of a type 
(usually the wild-type) that we call type 1. During therapy with drug A, the viral load decreases 
(eventually reaching undetectability) but in the long run a drug-resistant variant emerges, namely 
we observe the rise of the serum viral load of what we denote by type 2 virions. This forces us to 
change therapy and resort to drug B, which is able to strike type 2 down, but eventually another 



mutation arises (type 3 mutants), which escapes the action of drug B and serum viral load starts 
increasing again. We identify the moment when we detect type 3 as the moment of therapeutic 
failure (with the current tests at hand, this is usually when serum viral load reaches 20 copies/ml). 
Of course we may have at our disposal another drug capable of fighting this mutation too, and the 
process could continue, but we are interested in those patients which may eventually develop drug- 
resistance and a two-step mutation is the simplest case where one can test the difference between 
combination therapy and sequential monotherapy. Our main question is whether providing drug A 
and drug B together can significantly delay the instant of therapeutic failure (if so, it would be an 
argument in favour of combination therapy) . 

In Section [2] we describe a stochastic model (Definition I2.ip where, as in the basic model (jl.2p 
there is an interplay between free virus, infected hepatocytes and noninfected hepatocytes. The 
interactions are random and so are replications within the infected cells: with a small probability, 
whenever an infected cell produces a new virion, there might be a mutation, that is a random 
error in the replicative process. In particular when a cell is infected with the type 1 virus, it is 
possible that it also produces type 2 virions (more precisely we suppose that a mutation from type 
1 leads to type 2 virions), and, still with a small probability, type 3 is generated by mutation from 
hepatocytes infected with the type 2 variant. Thanks to the fact that the number of hepatocytes is 
very large, we are able to approximate our stochastic model with a deterministic model which, as 
the basic model, is the solution of an ODEs system, the system in (|2.4p . In Section [2] we recall the 
mathematical results which allow us to approximate the average behaviour of our stochastic model 
with this deterministic model (Definition 12.21 and Theorem 12. 3p and explain how to apply them to 
our specific case. 

In Section we analyze the equilibrium points of the ODEs system (|2.4[) (the system has no 
exact solutions). There are four equilibria: the disease- free equilibrium Xq, an equilibrium X\ 
with all the three types of infection, X2 with type 2 and type 3 infection together and a fourth 
equilibrium X3 with only the type 3 infection. The stability of these equilibria depends on the 
reproductive ratios of the three types (defined in (|3.5p ). In particular we are able to say that a 
variant tends to disappear if its reproductive ratio is smaller than 1 or if it is smaller than the 
reproductive ratio of another variant which has reproductive ratio larger than 1. We discuss how 
drug-resistance is explained by our model. More precisely, therapy lowers the reproductive ratios 
of the affected variants and thus, under therapy, the system moves from a situation X% to other 
equilibria where the cured variants disappear and others emerge. The main assumption is that, 
even if at X\ all three variants are present, nevertheless type 2 and type 3 are undetectable or at 
least they represent only a small fraction of the total viral load, hence they are not observed until 
the system reaches X2 or X3 (which are the "drug-resistant states"). 



In Section H] we focus on the main question of the paper, that is whether drug resistance appears 
first under combination therapy or under sequential monotherapy. The answer relies on the nu- 
merical solution of the ODEs system (I2.4p where therapy is considered in (I2.3P in the appropriate 
time intervals (see Remark \2A\i . To compute the numerical solutions we pick some of the unknown 
parameters from the literature and deduce the missing ones as functions of the serum viral loads at 
equilibrium (in Xi, X2 and A3) of the three variants. Letting these viral loads vary among several 
possible configurations (which depend on the fitness of the unknown variants type 2 and type 3) 
we study 46 plausible parameter sets and the corresponding numerical solution. Each solution 
is tested with combination therapy and sequential monotherapy. The time when type 3 reaches 
detectability under the two therapeutic approaches is computed. A comparison is presented, with 
different initial conditions. 

Section [5] is devoted to the discussion of the consequences of our results onto the choice of the best 
therapeutic approach in the cure of chronic HBV infection. In particular our results suggest that 
combination therapy is the best choice if started at the early stages of the chronic HBV infection. 

2. Stochastic model and deterministic approximation 

The simplest model for viral reproduction is the branching process (|22j). This model is very 
rough (the behaviour of the average number of particles is either exponential growth or exponential 
decrease), therefore to add complexity, one introduces space and/or interaction between particles, 
see |17| . Adding space one obtains the branching random walk where the reproductive capacity 
of the virus depends on its "location" (which may represent not only position but also type): its 
behaviour on complex networks has been studied for instance in [5l El [58]. In order to obtain 
a more realistic behaviour, where only a certain maximal viral load can be achieved (or at least, 
growth slows down when the viral population is too large) , we can put interaction into play (logistic 
competition is perhaps the simplest interaction). Spatially displaced interacting particle systems 
have been studied for instance in [21 [31 H] . 

Here we present a model, with no space, which is essentially a modification of the basic model 
(|1.2p . is stochastic and describes more complex interactions between the involved "particles". The 
actors in the basic model are T, Y and V, representing the total number of uninfected hepatocytes, 
of infected hepatocytes and of free virions, respectively. We first justify a modification of this model 
and introduce a stochastic description, then put different genomic variants into play. 

We assume that, by homeostasis, the total number of hepatocytes is constant, that is T-\-Y = Hq, 
where Hq is fixed. Namely we suppose that, whenever an infected hepatocyte dies, it is replaced 
(almost immediately) by an uninfected one. Of course this means that infected hepatocytes are not 
destroyed by the immune system at a too fast rate. This is in agreement with the common belief that 
cytotoxic T lymphocytes may directly inhibit viral replication and thus inactivate HBV without 
killing the infected hepatocyte, and that in chronically infected patients, infected hepatocytes escape 



immune recognition, sparing massive liver damage (see |12L H8|, [^Tj). This is particulary true during 
the first stage of the disease, the compensated phase (which lasts years under a correct antiviral 
therapy) or in the large group of HBV infected people called inactive carriers |47j . This assumption 
allows us to deal only with the dynamics of Y and V and to replace T with Hq — Y . Since the 
numbers Y and V change according to the random interactions between free virus, hepatocytes 
(infected or not) and immune system, it is natural to treat Y(t) and V(t) as random variables. 
More precisely, (Y(t), V(t))t>o is a continuous time Markov chain, with values in N 2 where N is the 
set of nonnegative integers (including zero). The random interactions are as follows 

(i) a single virion infects target cells at a rate which is proportional to the fraction of target 

hepatocytes among the total, namely a ^~ ; 
(ii) a single infected hepatocyte dies at rate d; 
(iii) virions are produced by any infected cell at rate p; 
(iv) virions are cleared from bloodstream at rate c. 




out of the system 



Figure 1. The dynamics of the one-type system. 
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5Y, if Y > 1; 

pY; 

cV, if V > 1; 

otherwise. 



Thus the dynamics of (Y(t),V(t))t>o is described by the following rates (q(k,l) denotes the tran- 
sition rate from k to I, k, I € N 2 ): 

'q((Y,V),(Y + l,V)) 
q((Y,V),(Y-l,V)) 
q((Y,V),(Y,V + l)) 
q((Y,V),(Y,V-l)) 
q(k,l) 

Now we introduce the three different variants of the HBV: the pre-existing variant type 1; type 
2 which appears after therapy with drug A as a mutant from type 1; type 3 which is a mutation 
from type 2, selected by therapy with drug B. We denote by V\ the total number of type 1 virions, 
V2 the total number of type 2 virions, V3 the total number of type 3 virions, and by Y\, I2, I3 
the number of hepatocytes that are carrying the corresponding type of infection (we suppose that 
an hepatocyte can be infected only by one virion at a time). The process is characterized by the 
following positive parameters: 



(a) the infection parameters a±, o<2, «3, where on corresponds to the rate at which virions of type 
% infect target hepatocytes (i = 1,2,3); 

(b) the production parameters p\, q2, r%, i.e. the rates at which hepatocytes which are infected 
with type 1, type 2 and type 3 virus respectively, produce virions of the same type; 

(c) the mutation parameters p2 and q%, i.e. the rates at which hepatocytes which are infected with 
type 1 and type 2 virus respectively, produce virions of type 2 and type 3; 

(d) the hepatocytes clearance rates S\ , 82 , £3 , where Si corresponds to the rate at which hepatocytes 
carrying type i infection are removed and thus replaced by target hepatocytes (i = 1,2,3); 

(e) the free virions clearance rates c\, C2, C3, where Cj corresponds to the rate at which virions of 
type i are cleared from bloodstream (i = 1,2,3). 

Definition 2.1 (Stochastic model). Given the positive parameters cti, 8i, Ci (i = 1,2,3,) and p±, 
Q2, ?"3; P2, Q3, the stochastic process (Yi(t),Y2(t),Y3(t),Vi(t),V2(t),Vs(t)) t >o is a continuous time 
Markov chain with values in N 6 . The transition rates qk,k+i (k € N 6 , I 6 Z, e ) are as follows, for all 
k^O, 
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otherwise; 



-e4, — es,— e6, respectively; 



(ei being the elements of the natural base of 

Although this is a random process, results by Kurtz prove that its dynamic behaviour is approx- 
imated by the solution of an ODEs system (|24|. [25], see also [15] for an example of application in 
epydemics models). These results apply to density dependent processes, hence we must consider a 
modification of our process. 

Definition 2.2 (Density dependent process and density process). A one-parameter family of con- 
tinuous time Markov chains (X^ N ' (t))t>o with state space E C TL d and transition rates (qij) is called 
density dependent if there exists a continuous function f : M. d x Z d — > R, such that 

q k ,k+l = N f(^ 1 )' l ^° and k > l e ^- 

Suppose (X^ N > (t))t>o is a density dependent process. By rescaling with N we get the density process 
(X N (t)) t > :=(±X( N )(t)) t > . 



Under certain conditions (-Xjv(£))t>o converges to a deterministic process that is the solution of 
a system of first order ODEs that is governed by the following function F: 

F(x) = Y / lf(x,l), 



as stated in [24^ Theorem 3.1]. We recall the theorem here. 

Theorem 2.3 (Deterministic Approximation). Suppose that there exists an open set E C M. d the 
function F is Lipschitz continuous and 

sup y \l\f{x, I) < oo; lim sup > \l\f(x,l)=0. 

i |(|>r 

Then, for every trajectory (x(s,xq),s > 0) satisfying the following system of ODEs 

—x(s,x ) = F(x(s,x )); 

x(0,xq) = 0, x(s,xo) £ E, < s < t, 
lim7v-i.cxD -Xjv(O) = xo implies that for every 5 > 0, 

lim P( sup \X N (s) - x(s,x )\ > 5) = 0. 

AT->oo 0<s<t 

This theorem implies that the process (-Xjv(£))*>o can be approximated to first order by a deter- 
ministic process, for large N. If the density process (X]y(t))t>o is initially close to xq, it will tend 
to stay close to the trajectory (x(s,xq),s > t) in some appropriate time-interval, subject to small 
random oscillations about the path. It is even possible to describe the behaviour of the random 
fluctuations of the density process (X7v(i))j>o around its deterministic approximation ([191 Ch.ll], 
see also [35]), but this goes beyond the aim of this paper. 

Now we apply these results to the stochastic model defined in Definition 12.11 If we take N = Hq, 

we have 

'«i(l -Xl -x 2 - x 3 )x i+3 

p 1 x 1 
f(x,l) = { p 2 x 1 + q 2 x 2 
Q3X 2 + r 3 x 3 

CX{ 

Thus F has the following expression and it is easily checked that the hypotheses of Theorem | 

are satisfied: 

(a\ (l-yi-y 2 - 3/3) vi ~ hyi\ 

"2 (1 - 2/1 - 2/2 - 2/3) V2 ~ 0~2V2 

«3 (1 - 2/1 - 2/2 - 2/3) ^3 - ^22/3 

PlVl ~ CVi 

P22/1 + 922/2 - CV 2 

\ 932/2 + r 3 y 3 - cv 3 ) 



if I = ei,i = 1, 2, 3; 

if I = —ei,i = 1, 2,3; 

if I = e^; 

if/ = e 5 ; 

if/ = e 6 ; 

if I = —ei,i = 4, 5,6; 

otherwise. 



F(yi,y 2 ,y3,vi,v 2 ,v 3 ) 



(2.3) 



Since Hq = 2 • 10 11 ([IS], also in accordance with the result of 139- 10 6 cells/g in [415] and an average 
weight for a human liver of 1.5-10 3 g), we may infer that the trajectories of the solution of the 
following system of ODEs 

(m,y2,m,vi,V2,v 3 ) T = F(y 1 ,y 2 ,y3,v 1 ,V2,v 3 ,). (2.4) 

are a good approximation of the ones of the stochastic process ( -^ , -jf-*- , -g-*- , -^ , -|M , -|M- J , 
at least in a time interval [0,£o] (where we assume that t$ is sufficiently large for our purposes). 
This means that y\{t), y2(t), ys(t) denote the fraction of hepatocytes that are infected, at time t, 
with type 1, type 2 and type 3 virus respectively. The total number of virions of type i (i = 1,2,3), 
flowing in the bloodstream at time t, is approximated by Vi(t) ■ 2 • 10 11 and the corresponding serum 
viral load is obtained dividing this quantity by 3 • 10 3 (the average number of ml of serum in the 
human body). 



Remark 2.4. The stochastic model we defined in (|2.ip and its approximation, that is the solution 
of (|2.4j) are not only models for the case where there is no therapy, but also for the ones where 
there is therapy with drug A or B or both. We only have to change the parameters involved by 
therapy. To be precise, in the first case we must replace a\, p\ andp2 by (1 — 1]a)c(i, (1 — £a)pi and 
(1 — ea)P2 respectively; in the second and third case a 2 , §2 and q 3 by (1 — t]b)oc2> (1 — £b)<12 and 
(1 — £b)q 3 respectively; in the third case a\, p\ and p2 must be replaced by (1 — r)c)ca, (1 — £c)pi 
and (l — £c)P2 respectively, where rjc and ec are the efficacy parameters of the combination of drug 
A and B against the type 1 infection. 

3. Analysis of the deterministic model 

The ODEs system (|2.4|) has no exact solution, thus in order to understand its behaviour, we 
start by analyzing its equilibria. Solving F(x) = (F taken from (|2.3p ). we find that there are 
four equilibrium points, the disease-free state and three infected equilibria. The Jacobian matrix 
of (|2.3p is J(x) given by 
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a 3 (3 


Pi 
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P2 


Q2 








-C2 





\ ° 


93 


r 3 








-c 3 ) 



where /3 = 1 - y x - y 2 - y$. 

Before going into the analysis of the equilibria, let us mention a few words about the choice 
of the parameters and what we expect to observe. The viral clearance parameters q depend 
only on the life-time of virions in bloodstream, thus on the antibodies production of the infected 
patient. We assume that this production does not depend on the variant, thus from now on we 
put c := c\ = C2 = c 3 . We also assume that P2 <C p\ and q 3 <C 02- Indeed P2 and ^3 represent 



the reproduction rate of type 1 and type 2 infected hepatocytes respectively, multiplied by the 
probability that a mutation appears (from type 1 to type 2 and from type 2 to type 3 respectively). 
The probability of mutation is very small, thus we assume that not only P2/P1 and (73/(72 are very 
small but also q 3 and P2/Q2 are small (in Section 0] we will choose (73 = 10~ 7 or (73 = 10 -6 and 
computation will lead us to, in different cases, P2/Q2 between 10~ 10 and 10 -4 ). Therefore in the 
analysis of some of the equilibria, we will approximate P2/Q2 and (73 with zero. An important role 
in the stability of the equilibria is played by the basic reproductive ratios of the three variants 

Ri = ^, Ba = ?f-, R * = ^ (3.5) 

C0\ C02 Cd 3 

Note that if the patient is undergoing therapy the corresponding infection and production param- 
eters have to be changed according to Remark 12.41 For instance, under therapy with drug A, the 
basic reproductive ratio of the type 1 variant becomes {a\p\{l — 77a)(1 — eaJ)/{c&x). 



3.1. Disease-free equilibrium. 


The disease-free equilibrium is 
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This matrix has six eigenvalues: 


-6 2 - c± 






-s 3 






—5\ — c ± s/iaipi + (5i - c) 2 


y / 4a 2 q2 + (82 ~ 


-c) 2 


- C ± y/4a 3 r 3 + (<*3 " 


-c) 2 



2 2 2 

These eigenvalues are real and they are all negative (hence Aq is locally asymptotically stable) if 
Ri < 1 for all i = 1,2,3. If at least one Ri is larger than 1, then there is at least one eigenvalue 
which is positive and Xq is locally asymptotically unstable. 

3.2. The infected equilibrium without type 1 and type 2 infections. The infected equilib- 
rium without type 1 infection is A3 = (0,0,y 3 ,0,0,v 3 ) and has biological meaning (that is, y 3 and 
v 3 are positive) only if R 3 > 1 , since 

£3 = (-R3 - 1)<V Q 3, 3/3 = 1 — I/R3 = a 3 v 3 /S 3 R 3 . 

The Jacobian matrix, evaluated at Ai, has six eigenvalues 



-a 3 r 3 (c 
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-a 3 r 3 - c 2 ± y / 4c 3 <5 3 + (a 3 r 3 - c 2 ) 5 



2c 



10 



They are all real and easy computations show that X$ is locally asymptotically stable if -R3 > 1, 
i?3 > R\ and R3 > i? 2 . If R3 < 1 or R2 > R3 or R\ > R3, then X3 is locally asymptotically 
unstable. 

3.3. The infected equilibrium without type 1 infection. The infected equilibrium without 
type 1 infection is X 2 = (0, y 2 , f/3, 0, v%, V3) and has biological meaning only if R2 > 1 and R2 > R3, 
since 

,Tn = 92 (fl 2 -l)(J?2-fl3) , 1q = 93(fl2-l) = 93 «2 .-;,_ 

1/2 c R2(R 2 -Rs+R3q3/r3) , 3 c(R 2 -i?3+«393/r 3 ) q2R2-R 3 2 

— c — — c /?■? — 

2/2 = ^2, 2/3 = FJ7?f V 3. 

Two eigenvalues are 

1 



- \-c - <5i ± Vpi ~ c ) 2 + 4aipi(l - y 2 - 2/3, 

which are negative if and only if R\ < 1/(1 — y 2 — 2/3) = R 2 - The other four eigenvalues are the 
solutions of a quartic equation and have expressions which are too involved to be studied here. 
We prefer to study the approximated solutions obtained substituting q% = (and 7)3 ~ 0, which is 
plausible if R2 — R3 is not too small) in the characteristic polynomial: 



- ( -a 2 v 2 - 62 - c ± y/(a 2 V2 + S 2 - c) 2 + 4a 2 g 2 (l - 2/2 - Vz, 



- (-c -8 3 ± y/(S 3 - c) 2 + 4a 3 r 3 (l - 2/2 - yz) 

Computations show that these approximate eigenvalues are all negative if R 2 > R3 and t> 2 > 0. 
Thus X2 is locally asymptotically stable if R 2 > 1 and R 2 > R3. If at least one of these conditions 
is not satisfied, X2 is locally asymptotically unstable. 

3.4. The infected equilibrium with type 1. The infected equilibrium where also type 1 is 
present, is X 1 = (#1, #2, #3, #1, £2,^3) where 

vi = s tvu v2 = Tt 2 m, h = ^m, 

f. _ H2(p2/g2) f . .-. _ P293 R2R3 V^-DJ 

« 2 ~ R1-R2 yi ' ^ 3 ~~ g 2 r 3 (R 1 ~R 2 )(Ri~Rs)- 

As for yi, it has a complicated expression which is of no interest here (but can be obtained from 
the relation y\ + y 2 + 2/3 = 1 — V-Ri)- Here is its approximation if one substitutes ^3 « 0: 

(i?i - l)(i?i - i? 2 ) 

2/1 ~ 



R 1 (R 1 -R 2 + R 2 ( P2 /q 2 )) ' 



thus Xi has biological meaning if i?i > 1, i?i > R 2 and i?i > i?3. The eigenvalues of J{X\) have 
too complicated expressions, hence we prefer to study the approximate solutions obtained with the 
substitutions qs = 0, p 2 = (and v 2 « 0, which is plausible if p 2 /</ 2 is small and Ri — R 2 is not too 
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small) in the characteristic polynomial. The approximated eigenvalues are 

- f-aiui -Sx-c± \/(aivi + 5[ - c) 2 + 4aipi(l - yi - y 2 - m) 

- (-c -5 2 ± y / (S 2 - c) 2 + Aa 2 q 2 (l - Vi 

- \-c - <5 3 ± y/(6 3 - 



2/2 ~ 2/3, 



+ 4a 3 r 3 (l -y\-y 2 - m, 



Substituting y\ + y 2 + y 3 = 1 — 1/R\, we get that X\ is locally asymptotically stable if R\ > 1, 
i?i > i?2 and Ri > i? 3 and unstable otherwise. 

The conditions for stability can be summarized by Table 1. 



Ri 


R2 


R3 


Equilibria 


Ri < 1 


R 2 <1 


^3 < 1 


Xq stable, other XiS unstable 


Ri > 1 


R 2 < R± 


i? 3 < i?i 


X\ stable, other XjS unstable 


R\ < R2 


R 2 >1 


i? 3 < i?2 


X 2 stable, other JQs unstable 


Ri <R 3 


R2<Rs 


R3>1 


X3 stable, other X^s unstable 



Table 1. Stability of the equilibria as a function of R±, R 2 , R3. 



3.5. Drug resistance explained by the model. Let us discuss how the presence of these equi- 
libria and their stability can be biologically interpreted. We are modelling the infection within a 
chronic patient, therefore with no treatment all three variants have reproductive ratios larger than 
1. Moreover we believe that type 1 has largest fitness and type 2 has larger fitness than type 3. 
Therefore in the beginning R\ > R 2 > R3 > 1. With no cure, after an appropriate amount of time, 
the system tends to lie in a neighbourhood of X\. We believe that the reason why drug resistant 
mutants are not observed unless under therapy is the fact that, even if mutations do happen, nev- 
ertheless in the equilibrium together with the wild type, mutants are numerically negligible. In 
other words, if we denote by Vij the i^-th component of the equilibrium point Xj, we assume that 
v 2 ,i and u 3> i are negligible if compared to vii (we already used this assumption in the analysis of 
equilibria). If the patient is given drug A, then type 1 is cured (R\ is lowered, by multiplication 
by (1 — tja) (1 — ea)), then the system moves towards equilibrium X 2 , where variant 2 emerges. If 
also type 2 is cured, with drug B (R 2 is lowered, by multiplication by (1 — ??b)(1 — £b)), the new 
equilibrium will be X 3 , that is, the patient will show drug resistance. 

Studies show that combination therapy has a larger efficacy in reducing the production rates 
(compared to single drug therapy); for instance [28] showed that, at the dosage of their trial, 
lamivudine only has e = 0.94, while combined with famciclovir it gets e = 0.988. We suppose that 
whenever a virus is the target of two drugs together, the effects add in an independent manner, 
namely the combined efficacies of the cocktail drug A + drug B against production of type 1 and 
against infection of new cells respectively, are Ec = 1 — (1 — £a)(1 — eb) = ea + eb — eaEb and 

r\c = 1 — (1— 77^4 ) ( 1 — tjb) = r/A+VB—VAVB- If the patient is given drug A and drug B in combination, 
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then the viral load of type 1 decays more rapidly (i?i is multiplied by (1 — r/c)(l — £q)) and the 
system moves directly from X\ to X%. The disease- free equilibrium represent the desirable state 
where the patient is permanently cured, and can be reached only if we have a drug at hand which 
acts also against the third variant (by assumptions, this is not the case under study). We believe 
that the fact that type 1 has the largest fitness is also reflected in the fact that when type 1 is 
removed from the system, the other two types can never exceed type l's initial viral load. Moreover, 
viral loads of a fixed type increase when the competitors are erased. Based on these observations, 
we make some more assumptions on the Vijs: v\,\ S> i>2,i > ^3,1, ^1,1 > ^2,2 > ^3,3 > ^3,2 > ^3,i> 

^2,2 > 1>2,1- 

The qualitative behaviour in time of the logarithm of viral loads under combination or sequential 
therapy (if therapy starts when the system is at X\) can be described by Figure 2 (where v\(t) 
is represented by a black line, V2(t) by a dashed line and vs(t) by the dotted one). Note that the 
picture is simplified, it does not show the biphasic decrease of viral loads under therapy, which is 
observed in trials and has been explained in |15j . 



via 




^2,1 
^3,3 



V3.1 




Figure 2. Viral loads dynamics under sequential (left) and combination (right) therapy. 



4. Numerical comparison under sequential or combination therapy 

The best therapeutic approach will be the one which sets off for a longer time the appearance 

of a detectable viral load of type 3 virions in the bloodstream. It is clear that the very presence of 

type 1 and type 2 virions in the bloodstream and in infected hepatocytes is an antagonist to the 

development of type 3, since all the variants are competing for a finite resource, namely, the target, 

non-infected hepatocytes. On the other hand, whenever type 1 (type 2 respectively) replicates, 

there is a chance of producing new variants (type 2 and type 3 respectively) that may be capable 

of drug resistance. Thus, in principle, it is not trivial to answer to the question whether it is a 

better strategy to employ sequential monotherapy (which may lead to surface type 2 virions and 

afterwards type 3 virions) or combination therapy (which strikes down both type 1 and type 2 but 

may clear the way to type 3). The answer relies on the quantitative dependence on time of viral 

loads. In order to analyze it, we resort to numerical solutions of the ODEs system (|2.4p . 
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We need to know the numerical values of the parameters involved, and of the initial condition. 
The twelve parameters are: the infection parameters a\, 02 and 03; the clearance parameters 
c, Si, 82 and £3; the production parameters pi, P2, q2, #3, T3. The initial condition is given by 
(yi(0), 2/2(0), 2/3 (0),ui(0),U2(0),U3(0)). Moreover, we also need the efficacy parameters rj and e of 
the drugs involved in the testing. We gather from the literature estimates on some of the parameters 
regarding the wild- type chronic infection and on the efficacy parameters. From [52J we have c = 
0.65/day, Si = 0.0143/day, e(adefovir) = 0.993, while from [33] we have e(lamivudine) ranging 
from 0.87 to 0.96 depending on the dosage. In their estimates, both papers use the assumption 
that a percentage varying between 5% and 40% of the total number of hepatocytes is productively 
infected in chronic HBV patients (the source of this assumption is [8]). Thus we take 2/1 (0) = 0.25, 
0.2 and 0.3 in different simulations. The order of magnitude of the initial viral load in the patients 
of [52] ranged from 10 7 to 10 9 copies/ml: for simplicity's sake, we take t>i(0) = 1 which corresponds 
to 2 ■ 10 11 free virions and a viral load of 0.67 • 10 8 (assuming that patients have 3 • 10 3 ml of serum). 

In our numerical tests, we choose ea = 0.97 and £b = 0.85: we assumed that the first drug 
is highly efficient against the wild type, while the second drug, which acts also on type 2, and is 
therefore more flexible, is less efficient (these efficacy parameters are nevertheless both in the range 
of the efficacies of known drugs). Although some papers assume that r/ is neither or 1 (see for 
instance [321 138] ), there are no explicit estimates of this parameter for the various drugs, hence in 
different simulations we pick tja and r]B in {0.25,0.5} ([32] used as a plausible value 0.5). Finally, 
it is not possible to gather information on the values of the infection and production parameters 
of the variants, since ideally we would like to compare sequential and combination therapy when 
facing any possible, unknown variants of the HBV virus. We assume that the clearance parameter 
of the infected hepatocytes does not depend on the type of virus which infected the cell, that is 
Si = 82 = 62, = 0.0143/day. 

In order to compute the missing parameters, we assume that we are given the numerical values of 
the Uj-components of the equilibria X\ (namely Via, 1^2,1 an d ^3,1), X2 (1^2,2 an d ^3,2)1 an d X$ (1^3, 3). 
Since these six equilibrium viral loads depend on the six missing parameters (see the expressions 
in in Section [3]), we take the inverse function and derive a±, 02, 03, P2, 12 and r%. The parameter 
p\ is derived from pi = M , once we fix the value of 2/1,1 (i/ij being the 2/j-component of X,). We 
fix 2/1,1 = 2/1 (0)- We are left with an arbitrary choice of q^, so we assume that it is very small and 
choose it equal to 10~ 7 . 

We proceed computing all the parameter sets which are obtained when the equilibrium viral loads 
vary among plausible values (we test for each value different powers of 10). We do not consider 
values that lead to nonpositive parameters (among the conditions there is ^3,2 > t>2,2^-); moreover 
we require that pi > q2 > ^3 (type 1 has a faster production rate, and type 2 is faster than type 3). 
The resulting ranges, which reflect our hypotheses on the equilibria (for instance the negligibility 
of 1*2,1 and 1*3,1), are listed in Table 2. We also require that that, at equilibria, the percentage 
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«2,1 


U3,l 


^2,2 


«3,2 


^3,3 


min 


io- y 


10- 9 


lO"" 


10^2^ 


io _b 


max 


min (10^, ^f) 


U2,l 


1 


/ '«2,2 "3,3 \ 


1 



Table 2. Range of equilibrium viral loads. 



of infected hepatocytes lies between 5% and 40% (parameter sets which do not satisfy all these 
conditions are disregarded). Under these conditions we get 46 cases to analyze, which correspond 
to different possible type 2 and type 3 variants (with more or less ability to infect and reproduce). 

Once we obtain the parameters, we numerically compute the solutions of the ODEs system in 
case of combination therapy (drug 1 and drug 2 given together until type 1 and type 2 are both 
eradicated) or of sequential monotherapy (drug 1 is given until type 1 is cured and we start giving 
drug 2 when type 2 reaches a viral load of 667 copies/ml). Combination therapy is judged as the 
best therapeutic approach if the time when viral load of type 3 exceeds 20 copies/ml is significantly 
larger than the same time under sequential monotherapy. 

In the numerical computations the arbitrary parameters such as r/A, tjb, yi(0) and q% have been 
varied (one at a time) with no significant changes in our conclusions. The only values on which the 
best therapeutic approach depends are the initial viral loads of type 3. We assume that the initial 
viral load ^(O) is 7 or 67 copies/ml and 2/2(0) = ^^(0) (that is, the fraction of initially infected 
hepatocytes carrying the type 2 infection is the same function of the corresponding viral load as 
j/2,1 is a function of t>2,i)- 

We sum up the conclusions of our numerical simulations in Table 3. The answer to our main 

question depends on the values of 1^3(0) and 2/3(0), which are the initial serum viral load and the 

fraction of infected hepatocytes of type 3 respectively. In the second column (CT best) we put 

the fraction of cases where the combination therapy is advantageous, and among these cases we 

compute the minimum, maximum and average percentage delay (mind, maxd and d respectively). 

Finally, T s and T c are the average times before drug resistance (i.e. the detection of type 3) for 

patients under sequential or combination therapy respectively. The table refers to data relative to 

the 46 cases obtained with the following choice of the arbitrary parameters: q% = 10~ 7 , 2/1 (0) = 0.25, 

t>2(0) = 10~ 7 , tja = i]B = 0.5. We changed the value of these parameters, one at a time, without 

significant changes in the conclusions, which is an indication of robustness of our comparison. It 

appears that the advantage of combination therapy is appalling if started before that type 3 is 

present (first row of Table 3). If therapy starts when type 3 has a viral load of approximately 1 

copy/ml but has not infected any hepatocyte yet (second row of Table 3), combination therapy is 

still the best choice in 38 cases among 46, but delayed drug resistance by more than 10% of the 

time only in 18 cases. Finally, if infection of type 3 has already progressed we look at the third row 

of Table 3. There we chose 2/3(0) = -^§-^3(0) (the fraction of infected hepatocytes corresponding 
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to type 3 in equilibrium together with the other strains, if ^3(0) coincided with V3, see (|3.6|) ) and 
obtained that while combination therapy is preferable in 16 cases out of 46, sequential monotherapy 
is preferable in 9 cases out of 46. In the remaining 21 cases there is substantial parity of outcomes 
under the two therapeutic approaches. 



vl 3 (0), y 3 (0) 


CT best 


minti 


maxd 


d 


f s 


f c 


0, 


46/46 


33% 


90% 


50% 


6.3 years 


9.5 years 


1,0 


38/46 


1% 


56% 


15% 


5.1 years 


5.8 years 


1, >0 


16/46 


0.4% 


7.7% 


3% 


2.87 years 


2.91 years 



Table 3. Comparison of therapeutic approaches depending on the initial conditions. 



5. Discussion 

We developed a stochastic model of HBV viral dynamics (with or without therapy) which ac- 
counts not only for different aspects of the interplay between the virus and the host, through some 
parameters, but also for the presence of heterogeneous strains in the viral progeny, due to random 
errors in the replicative process. In particular we considered two consecutive mutations of the 
primary infection. The numerical value of the parameters in the model depend on the patient 
and on the viral strains. Indeed the host is responsible of an immune response that counteracts 
the infection through viral clearance from the bloodstream (humoral immunity, parameter c) and 
clearance of infected hepatocytes (cell-mediated immunity, parameters 5i). On the other hand, 
any virus has a specific capacity of infection of target cells (parameters a,) and of reproduction 
through the infected cells (parameters p\, q2 and r% for the three variants, respectively). More- 
over whenever there is production of type 1 virus there is a small probability of producing a type 
2 virion and analogously, from type 2 infected hepatocytes there might be production of type 3 
virions (the mutation rates are P2 and q% respectively). We approximated the behaviour of the 
stochastic model with the behaviour of a deterministic model, given by the solution of the ordinary 
differential equations system (|2.4p . which was then the object of our study. 

Under the natural assumption R± > R2 > R3, based on the fact that more easily observed 
variants must have a larger fitness (reflected on the order of reproductive ratios) , the model is able 
to explain drug resistance. Indeed, as we showed in Section [3l if there is no therapy then in the long 
run the mutants, even if they appear by mutation, remain numerically negligible. If the patient 
undergoes therapy which removes only the primary virus, then the drug-resistant strain type 2 
surfaces; while if therapy acts also against this latter strain, in the long run the third variant, 
which by definition we suppose drug-resistant to any drug at hand, appears. Note that the model 
gives a qualitative description of the appearance of drug-resistant strains even if the characteristics 
of these strains are unknown (the only assumption is that mutants are negligible when together 
with the wild- type and in general have smaller fitness). 
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We used this model to compare the efficacy of combination therapy versus sequential monother- 
apy. In this comparison, the best therapy is the one that mostly delays the appearance of the drug 
resistant variant (i.e. the time when this variant reaches a serum viral load of 20 cp/ml). The 
analysis of the equilibria (Section [3]) tells us how the system behaves qualitatively with or without 
therapy but does not tell us how long it takes to reach any of these equilibria. This can be deter- 
mined only solving f|2.4[) , where the infection and production parameters are multiplied by 1 — fj and 
by 1 — e respectively, and f\ and e depend on the therapy. For instance, in sequential monotherapy 
for t € [0,£i] therapy is drug A, for t > t\ therapy is drug B, where t% is the time where type 2 
reaches a worrying viral load (£i depends on the solution itself). The system can be solved only 
numerically (through a mathematics software, we used Maple). To plug in numerical values of the 
parameters, we took c±, 5\ from the literature and assumed c\ = C2 = C3 and 6\ = 82 = 83. Six 
unknown parameters (qi, 0:2, 03, P2, Q2, ^3) were obtained from the expression of the equilibria, 
as functions of the viral loads at equilibrium of the three types. The parameter p\ is a function 
of the equilibrium viral load of the primary infection and of the corresponding fraction of infected 
hepatocytes (which were gathered from the literature); for the last unknown parameter, (73, we had 
to make an arbitrary choice (and put it equal to 10~ 6 or 10™ 7 in different simulations). We believe 
that, as long as qs is small, its numerical value can only rescale the time of emergence of the third 
variant (the smaller q^, the longer the time) but should not change the kind of therapy which is 
most suitable. As for the choice of the efficacy parameters of the drugs involved in the testing, we 
picked them from the literature. 

We note here that under the assumption that 5± = 62 = 63, for all plausible choices of the 
equilibrium viral loads, we obtained a\ < 02 < 03. This means that, if the mutants are not 
more capable of the wild-type to evade the cell-mediated immune response, then only a greater 
ability of infecting target hepatocytes can lead them to detectability. We performed some numerical 
calculations with some choices of 82 and £3 such that 8\ > 82 > #3 and obtained that in these cases 
it is possible that a% > 02 > «3 (which is more intuitive, since one is lead to believe that larger 
fitness also means greater infection ability). Thus, in our model, drug-resistant mutants can reach 
detectibility if, compared to the wild-type, they either have a greater infection ability or if they are 
better at masking themselves into infected cells. This is plausible, since mutations may for instance 
affect epitopes recognized by cell-mediated immune response. 

It is worth noting that, even if we were forced to make some arbitrary parameter choices and 
the model is not fitted to any data (which is intrinsic in the problem itself since we want to address 
unknown variants of the HBV), nevertheless the order of magnitude of the time before drug- 
resistance that we obtained is in accordance with the time observed in case trials. For instance in 
|29l Table 1] many patients report, with different drugs, resistance after a treatment period ranging 
from 3 to 6 years. 
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Our numerical analysis shows that combination therapy is by far the best choice if started at the 
early stages of the chronic infection, while it seems that the advantage is not striking if therapy is 
in act when drug-resistant mutants are already numerous and have started to significantly infect 
the hepatocytes. This difference is not surprising, since if type 3 mutants are absent, they are 
generated only by mutation from type 2, hence the need of lowering the number of type 2 virions 
and infected cells. On the other hand, if type 3 mutants are already present, in some cases the 
competition with other strains might slow the proliferation of the type 3 infection. Nevertheless, we 
believe that in reality early therapy is more convenient also for other reasons. In our deterministic 
approximation of the model, there is a constant production of mutants, in a fixed proportion of 
the existing replications. In reality those replications are random and the time r of appearance of 
mutants is random. If there are many replications r is more likely to be small; therefore a good 
strategy is also to keep replications low in order to delay the first appearance of mutants. This 
randomization is not considered in the deterministic approximation of the model and gives another 
indication in favour of early combination therapy. The biological and clinical importance of this 
model is to show a clear advantage in favor of the combination strategy, as it could lead to a greater 
delay in the development of drug resistant variants, which have less fitness than the wild type (drug 
sensible) strain, but are able eventually to escape from the therapy. From the analysis of the model 
we also want put emphasis on the need of an early antiviral therapy: since the probability of the 
emergence of mutation is directly proportional to the viral fitness, high levels of wild strain virus 
could generate more easily drug resistant variants (type 3 mutants). Consequently the probability 
of mutation, over time, is inversely proportional to the selective pressure exerted by the antiviral 
therapy. 

By the way, there is also clinical evidence that the best choice is a broad spectrum early therapy: 
for instance in [50] the authors observed that the emergence of entecavir-resistant strains is more 
likely to appear in patients which were already lamivudine-resistant (rather than in treatment-naive 
patients). 

Concerning the type of combination therapy, ideally the best combination therapy for HBV 
infection should consider the association between a high power drug (with a laevorotatory structure) 
and a high genetic barrier drug (with a flexible, plastic structure). Power is the speed with which 
a drug causes the suppression of viral replication hence it is proportional to the efficacy parameter 
e introduced in Section H|). We can assume that a potent drug is a molecule that has a high affinity 
to the target protein. Such affinity, in addiction to an intrinsic chemical characteristic, could be an 
expression of a better adaptation to the chirality of nature of the target. Respecting the symmetry 
of nature, the laevogyre molecules are more selective and more powerful to the target, with fewer 
number of side effects. The genetic barrier represents the number of mutations needed by the virus 
to replicate effectively in the presence of the drug (a higher genetic barrier indicates that there 
is a smaller set of mutants that evade the action of the drug). Since the action of a drug is an 
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expression of the interaction active molecule - target protein, we can assume that a drug with an 
high genetic barrier is a molecule that has a "high flexibility" , which enables to adapt to different 
conformations of the target protein induced by punctiform genotypic mutations. That synergism 
is something more than the simple sum of the effects of two drugs, and it contains a stronger 
action than monotherapy, even with entecavir or tenofovir, molecules that to date represent the 
best compromise between power and genetic barrier. These two drugs show a low incidence of 
mutations after 4 years of treatment: 1.2% and 0% respectively, [29j (although in principle sooner 
or later drug-resistant mutants may appear with any drug, which is one argument in favour of 
combination therapy even with these drugs at hand). As we already mentioned this dramatically 
changes with entecavir in non-treatment-naive patients ([50j). Moreover recent results show the 
appearance of drug-resistance in coinfected HBV-HIV patients, treated with tenofovir ([36]). The 
authors identify in poor compliance to therapy the reason of this failure; it is our opinion that this 
is another aspect to keep in mind. In conclusion, basing on the numerical evaluations carried out 
in the model, we suggest to intervene in the natural history of HBV infection immediately and 
with a broad-spectrum approach: a combination of powerful and high genetic barrier drugs. We 
believe that combination therapy is less likely to select resistant strains, especially in patients with 
a suboptimal compliance. We may infer the same conclusions also considering the evidence-based 
effectiveness of the combination therapy (compared to the monotherapy) in other viral infections 
supported by multiple subtypes, termed viral quasispecies, such as hepatitis C virus (HCV) and 
human immunodeficiency virus (HIV) infections. 
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